Paternal hypoxia exposure primes offspring for increased hypoxia resistance

Background In a time of rapid environmental change, understanding how the challenges experienced by one generation can influence the fitness of future generations is critically needed. Using tolerance assays and transcriptomic and methylome approaches, we use zebrafish as a model to investigate cross-generational acclimation to hypoxia. Results We show that short-term paternal exposure to hypoxia endows offspring with greater tolerance to acute hypoxia. We detected two hemoglobin genes that are significantly upregulated by more than 6-fold in the offspring of hypoxia exposed males. Moreover, the offspring which maintained equilibrium the longest showed greatest upregulation in hemoglobin expression. We did not detect differential methylation at any of the differentially expressed genes, suggesting that other epigenetic mechanisms are responsible for alterations in gene expression. Conclusions Overall, our findings suggest that an epigenetic memory of past hypoxia exposure is maintained and that this environmentally induced information is transferred to subsequent generations, pre-acclimating progeny to cope with hypoxic conditions. Supplementary Information The online version contains supplementary material available at 10.1186/s12915-022-01389-x.

Hypoxia, defined as sufficiently low levels of oxygen to deprive tissues of oxygen, is a major physiological challenge [39][40][41], precipitating conserved physiological effects in a wide array of vertebrates [42][43][44]. The aerobic lifestyle of most animals requires a constant supply of sufficient oxygen, and low oxygen levels constitute a major environmental threat [45]. Hypoxia is particularly wellstudied in aquatic species, where physiological stress of hypoxia challenges many freshwater and marine organisms, and hypoxia is one of the most widespread issues in aquatic habitats due to the rise of dead zones and climate change [39,46]. Thus, it is now more important than ever to understand how different aquatic species will respond to environmental hypoxia.
The negative impacts of hypoxia on growth and reproduction are well-recognized [44,47,48], with recent studies also discovering transgenerational effects, whereby marine medaka fish (Oryzias melastigma) exposed to hypoxia show reproductive impairments in F1 and F2 generations, suggesting that hypoxia might pose a long-lasting threat to fish populations [40,49]. Fishes are notable for their adaptive abilities to acclimate to hypoxic conditions [39,[50][51][52][53], including physiological, morphological, and behavioral phenotypic responses. Remarkably, goldfish (Carassius auratus) and crucian carp (Carassius carassius) exhibit gill remodeling in response to hypoxia exposure [54,55]. Zebrafish (Danio rerio) embryos and larvae modify cardiac activity and blood vessel formation in response to hypoxic exposures [56]. In addition to morphological alternations, gene expression can also be affected by hypoxic conditions. In zebrafish gill tissue, more than 300 genes are differentially expressed between hypoxia exposed individuals versus controls; these changes in gene expression are correlated with morphological changes in gill structure, such as increased surface area of gill tissue [57]. Prolonged exposure to low oxygen has been shown to improve hypoxia tolerance in Murray cod (Maccullochella peelii [58] and snapper (Pagrus auratus [59]). Likewise, pre-acclimation of zebrafish larvae to mild hypoxia significantly improves their resistance to lethal hypoxia, and upregulation of some oxygen transport genes are associated with this acclimation [41]. Intriguingly, zebrafish offspring of males and females exposed to chronic hypoxia have a higher resistance to acute hypoxia than those of controls [25], despite the offspring having never been exposed to hypoxia, suggesting that parents may pass on information that may pre-acclimate offspring to cope with hypoxic conditions.
Gene expression and transgenerational effects can be regulated by epigenetic modifications, including DNA methylation, histone modifications, and non-coding microRNAs [3]. For example, in the marine medaka studies, transgenerational effects are associated with differential methylation in sperm and ovary, with altered gene expression in genes known to be associated with spermatogenesis and gene silencing [40] and cell cycle control and cell apoptosis [49]. DNA methylation has also been investigated as a means underlying transgenerational effects of chemical exposures in zebrafish [60], and the paternal methylome is believed to be stably transmitted to offspring in zebrafish [61,62] without global reprogramming in primoridal germ cells [63,64], potentially facilitating environmental specific information transfer.
Here we use zebrafish (D. rerio) to further explore the phenomenon of transgenerational acclimation of hypoxia in fishes. We focus on paternal exposure as we predict that environmental specific information can be transferred via sperm, as observed in other studies [4,40,65,66]. We test whether paternal exposure to hypoxia stimulates phenotypic responses in offspring, using behavioral phenotyping to identify resistance to acute hypoxia (time to loss of equilibrium). We then use RNA-Seq to identify candidate genes that are differentially expressed in control and hypoxic progeny and whole genome bisulfite sequencing data to assess whether changes in DNA methylation underpin alterations in phenotype and gene expression (Fig. 1).

Hypoxia tolerance assays in unexposed offspring
Progeny of males exposed to moderate hypoxia for 14 days show a greater resistance to acute hypoxia than progeny of control males. Time to loss of equilibrium was, on average, 32 s longer for progeny of hypoxia exposed males (t = 2.52, p = 0.014; confidence interval (CI) = 7.38, 58.83; Fig. 2A) and the progeny of hypoxia exposed males lost equilibrium, on average, 9.18 times vs. 12.32 times in the control progeny (t = − 4.21, p <0.0001, CI = − 0.433, − 0.158; Fig. 2B). However, when accounting for family (to avoid pseudo-replication due to multiple offspring per family being tested), the effect of time to loss of equilibrium becomes nonsignificant (t = 1.67, p = 0.133, CI = − 12.66, 79.22; note this model accounts for observed heterogeneity in the data), though the effect remains for the number of times offspring lost equilibrium (z = − 2.46, p = 0.014, CI = − 0.567, − 0.037). Hence, there are strong family effects (time to loss of equilibrium Treatment ID variance = 42.38; loss of equilibrium frequency Treatment ID variance = 0.025), with greater resistance to acute hypoxia being observed for just two families, H1 and H2 (Fig. 2C, D).

Transcriptome wide gene expression patterns in unexposed offspring
We sequenced 17.4-20.0 million reads per sample. Mapping assigned 14.1-16.1 million reads to genes, representing 80.15-81.41% of reads uniquely assigned to a gene, with detected expression in 26,260 of the 32,057 genes in the GRCz11.99 annotation.
A total of 91 genes were significantly differentially expressed between the offspring of control and hypoxia exposed males ( Fig. 3; Additional file 1: Table S1). Of the 91 differentially expressed genes, 36 genes were significantly upregulated, and 55 were significantly downregulated in the offspring of hypoxia exposed males (Fig. 3). Eight genes exhibited greater than 4-fold change in differential expression (2 upregulated in hypoxia, 6 downregulated in hypoxia) (Fig. 4). Most notably, two hemoglobin genes (hbaa1 and hbz (si:ch211-5k11.6); Figs. 4 and 5A) were upregulated by more than 6-fold in offspring of hypoxia exposed F0 males. Three additional genes were upregulated by more than 4-fold, but did not pass false discovery rate (FDR) correction-these included a third hemoglobin gene (hbba1; Figs. 4 and 5A; Additional file 1: Table S1) and a major histocompatibility gene (mhc1ula; Fig. 4). Importantly, the observed upregulation in hemoglobin gene expression correlates with the differential family effects observed in our loss of equilibrium tolerance assays (compare  Table S2). Sequenced offspring from the H1 and H2 families, the two families that show the greatest tolerance to acute hypoxia, also shows the greatest upregulation in gene expression. Indeed, the sequenced H1 offspring never lost equilibrium within our 240s long assay and the sequenced H2 offspring first lost equilibrium at 216.5s, whereas the third hypoxia offspring, from family H3 (with a first loss of equilibrium at 181.5s), and the offspring from C1-C3 (first loss of equilibrium > 3 s at 38.8, 24.2 and 18.9s, respectively) show little to no expression at these hemoglobin genes.
Six genes were downregulated by more than 4-fold in offspring of hypoxia exposed males (Fig. 4), including timm23b (translocase of inner mitochondrial membrane 23 homolog b (yeast)), acot21 (acyl-CoA thioesterase 21), znf1156 (zinc finger protein 1156), cyt1l (type I cytokeratin), irge1 (immunity-related GTP-ase), and zgc:158417. None of these downregulated genes appear to associate with hypoxia in any way, and hence, are not considered further. Five males from each treatment were then used to create F1 progeny, crossing the males to unexposed females (n = 5 families/treatment), with half of the sperm used for whole genome bisulfite sequencing, to assess differential methylation (n = 3 per treatment). At 20-21 days post fertilization, offspring undergo acute hypoxia tolerance (6-8 offspring/family) assays and n = 3 offspring/treatment are used for differential gene expression analysis (RNA-Seq) Of the genes we analyzed, 13,698 were associated with a GO term (including 58 of the 91 differentially expressed genes). Nine GO terms were overrepresented in our differentially expressed genes at a q-value cutoff of 0.1 (Table 1). Several significantly enriched GO terms were associated with lytic or proteolytic activity (i.e., serine hydrolase activity, serine-type endopeptidase activity, hydrolase activity, serine-type peptidase activity, glutamate decarboxylase activity, endopeptidase activity, peptidase activity, peptidase activity, acting on L-amino acid peptides). Proteolysis (GO:0006508 p-value = 7.43 × 10 −04 ; q-value = 0.7640) tended to be over-represented. Genes associated with the process of aging (GO:0007568; p-value = 1.99 × 10 −05 ; q-value = 0.1640) as well as to the response to oxidative stress (GO:0006979; p-value = 9.20 × 10 −04 ; q-value = 0.841) also tended to be over-represented.

The relationship between methylation in parental sperm and gene expression in offspring
To explore potential mechanistic explanations underlying increased tolerance to acute hypoxia and differential gene expression in the zebrafish offspring, we conducted whole genome bisulfite sequencing (WGBS) of sperm from hypoxic (n = 3) and control males (n = 3). We obtained an average of 67,891,334 reads per methylome (SEM ± 2,551,964) and a mapping efficiency of 45.02% (SEM ± 0.07%). Global cytosine-guanine (CG) dinucleotide methylation levels showed no differences between hypoxia and control samples (84.15% and 84.23%, respectively; p > 0.05; Additional file 3: Table S3). . Boxes illustrate the interquartile range (medians, 25th and 75th percentiles), and whiskers illustrate 1.5 * the interquartile range, above and below the 75th and 25th percentiles. Statistical analysis: A p = 0.014, B p < 0.0001 treatment effects modeled using generalized linear models. C p = 0.133, D p = 0.014 treatment effects modeled using generalized linear mixed models with family ID included as a random effect  In jawed vertebrates, methylation of transcription start sites (TSS), and in particular regions enriched on CG dinucleotides (CGIs), is associated with transcriptional silencing [67]. For two example cases, the upregulated hemoglobin cluster and the downregulated acot21 gene, we characterized methylation at individual positions and CGI using statistical and experimental criteria [68,69]. No differences in methylation levels were observed for these genes (Fig. 5A). Next, we explored methylation levels for all differentially expressed genes. Using a threshold of 20 methylation calls, we obtained 29 and 23 genes for the upregulated and downregulated groups. Whereas downregulated genes are hypermethylated and upregulated genes are hypomethylated for hypoxia and control samples, TSS methylation values remain stable when both groups of samples are compared (Fig. 5B). Finally, we explored the global effect of DNA methylation on gene expression using expression quantiles. We found methylation levels in the parental sperm and gene expression in the F1 offspring are coupled in hypoxia and control samples (Fig. 5C). Nevertheless, when we completed a differential methylation analysis of all genes, we found 420 significantly differentially methylation regions (226 hypermethylated and 194 hypomethylated) in response to low-oxygen treatment (Additional file 4: Table S4); however, these were not associated with corresponding alterations in gene expression in the offspring.

Discussion
We demonstrate that paternal exposure to hypoxia alters both the phenotypic response to hypoxia and gene expression in the offspring. Larvae of fathers that experienced moderate hypoxia-maintained equilibrium in acute hypoxia for longer than those of controls, indicating that paternal exposure stimulated a higher tolerance to hypoxic conditions. Using next-generation sequencing, we also detected significant changes in gene expression between control and hypoxia offspring, with two key hemoglobin genes upregulated in offspring of hypoxia exposed males-genes which may mediate the observed phenotypic differences, as they are involved in oxygen transport. This pattern of inheritance, through the paternal line, could have large evolutionary consequences as fathers are able to pass down valuable information to offspring that may enable better survival. However, the underlying mechanism for this transmission remains unknown, as we did not detect any differential methylation in the sperm of parental males at any of the differentially expressed genes. Importantly, our study provides evidence of increased tolerance to acute hypoxia through paternal exposure, thus, adding to increasing evidence that environmental challenges experienced by ancestors can provide progeny with environmental specific information that might allow future generations to survive the same environmental challenge, i.e., transgenerational plasticity [16,17,23,27,[29][30][31]33], though see [70] for an amphipod study that did not find evidence for transgenerational plasticity to hypoxia.
In a previous zebrafish study, 20 dpf larvae exposed to 4 kPa pO 2 revealed that offspring of parents exposed to hypoxia had longer time to loss of equilibrium (hypoxia resistance [25];). In preliminary trials, we found no behavioral differences using these oxygen parameters; thus, we increased the magnitude of hypoxia, exposing larvae to ~1 kPa, which produced the expected phenotypic differences. Importantly, the previous study [25] did not differentiate between maternal and paternal effects and did not provide a possible mechanism underlying this phenotypic effect.
It is important to note that our study detects intergenerational acclimation, with potential for transgenerational acclimation. A transgenerational study requires rearing fish through to create an F2 generation. Indeed, few studies that claim transgenerational acclimation are truly transgenerational, as the studies are conducted across a single generation [3]. Even if parents are exposed to environmental challenges prior to maturity, the primordial germ cells are still exposed to the challenge as well. Regardless, our results suggest that parents are passing on environmental specific information that may benefit offspring survival, thus facilitating acclimation to environmental conditions. Additional research is needed to investigate whether the offspring of hypoxia exposed parents also exhibit any changes in gill morphology, metabolic or cardiac function or behavioral responses, like surface respiration, to hypoxia, as observed in adult fishes exposed to hypoxia [51]. To try and understand the underlying mechanism priming progeny to better cope with hypoxic conditions, we conducted transcriptomic analysis of the progeny.
We detected 91 genes that were differentially expressed in the offspring of paternal males that were exposed to moderate hypoxia for 2 weeks. Most notably, two hemoglobin genes (hbaa1 and hbz) exhibited over 6-fold differential expression, were upregulated in offspring of males exposed to hypoxia, with another hemoglobin gene, hbba1, also upregulated, though non-significantly. Remarkably, sequenced offspring from the H1 and H2 families, the two families that show the greatest tolerance to acute hypoxia also show the greatest upregulation in hemoglobin gene expression (Additional file 2: Table S2). hbaa1, hbz, and hbba1 are all found on chromosome 3, are part of the major hemoglobin enhancer locus responsible for heme binding, and are instrumental in oxygen transport [71]. Exposing fish to hypoxia is typically considered to improve hypoxia tolerance through alterations in hemoglobin, hemoglobin-O 2 binding affinities, or cardiac function to improve low O 2 performance [59]. Thus, differential expression of hemoglobin genes (upregulated in offspring from hypoxia exposed fathers) could indicate altered physiological mechanisms to combat low oxygen, which could be precipitating the increased tolerance to acute hypoxia in the H1 and H2 families.
The relationship between hypoxia and hemoglobin function has been studied in larval zebrafish, suggesting that zebrafish larvae might be able to upregulate hemoglobin concentration in response to chronic hypoxia [72]. Under normoxic conditions, oxygen supply via diffusion seems to be sufficient to meet metabolic demands up to 12-14 dpf in zebrafish, but larvae appear to be able to use a circulatory system as a backup where necessary [72]. Further, impaired hemoglobin function does not impair routine oxygen usage in normoxia or at moderate levels of hypoxia in 5-42 dpf larvae, but functional hemoglobin does allow larvae to sequester extra oxygen from water in extreme hypoxic conditions [73]. The upregulated hemoglobin genes identified in this study are adult hemoglobins-larvae express embryonic/larval globins early in development, but somewhere between days 16 and 22, the embryonic/larval globins begin to decline [71,74], and adult globin expression increases, with the adult globin expression pattern nearly completely established by day 32dpf. hbae5 (ENSDARG00000045142) is the only embryonic/larval hemoglobin to show any pattern of differential expression in our transcriptomic data. We found that this gene was highly expressed and upregulated twice as much in our hypoxia offspring than in our control, though not significantly differentiated after FDR correction (p = 0.014, q = 0.391). hbae5 expression appears to peak around 22 dpf [71], so upregulation of this gene should result in higher affinity for oxygen, which may help explain the increased tolerance to the acute hypoxia that we observed in our 20-21 dpf progeny. hbae5 has also been shown to be upregulated more than 5-fold by hypoxia in zebrafish larvae directly exposed to hypoxia; ( [41]; note that hbae5 is called hbz in their study).
In addition to two genes upregulated over 6-fold in offspring of hypoxia exposed males, six genes were downregulated by 4-fold as well. These include timm23b (translocase of inner mitochondrial membrane 23 homolog b (yeast)), an integral component of membranes. Gene acot21 (acyl-CoA thioesterase 12), is a key component of acyl-CoA metabolic processes with thiolester hydrolase activity and found in the cytoplasm. Gene znf1156 (zinc finger protein 1156) has metal ion binding functions. Gene zgc:158417 and irge1 (immunityrelated GTP-ase family) are predicted to be integral component of membranes and used in GTP binding. Gene cyt1l (type I cytokeratin) is predicted to have structural molecule activity. Little information on these genes is available, and it is not clear how they link with hypoxia in anyway. Expression changes of such significant magnitude are likely to have large consequences, so some additional research around hypoxia and gene function is needed.
An HRGFish database [75] reports 50 key genes that are altered through hypoxia, but we did not find any of them to be differentially expressed in our transcriptomic data (Additional file 1: Table S1, showing data for 46 of the 50 HRGFish genes referenced for the zebrafish), suggesting that the genes that are altered by direct exposure to hypoxia may be different to those that may be passed on to future generations of hypoxia exposed ancestors. Importantly, transcriptomic studies of F1 and F2 adult tissues, like heart, gill, brain and liver, might reveal differential expression of some of these key hypoxia genes. Interestingly, our GO analysis highlighted an overrepresentation of genes associated with aging and oxidative stress. Studies have previously hypothesized a link between hypoxia (and effects on respiration), oxidative stress (free radical production) and aging in humans, with several studies demonstrating age-related changes in the hypoxia inducible factor system [76,77].
DNA methylation is known to be altered by hypoxia exposure in fishes [40,49] and other vertebrates [78,79]. Thus, we hypothesized that hypoxia might produce changes in the paternal methylome, which are inherited to the offspring, explaining the observed tolerance to acute hypoxia and differential gene expression. Interesting, we found methylation levels and gene expression were coupled in our hypoxia and control samples, despite the different tissues of origin for the methylome and RNA expression data. This coupling supports a common pattern between the paternal methylome and expression levels of the embryo. Yet, we observed no differences in methylation levels in any differentially expressed genes, including the hemoglobin cluster, suggesting that sperm methylation patterns are not responsible for changes in the observed changes in gene expression in our study, or perhaps through epigenetic control away from the focal genes. The medaka studies [40,49] exposed fish to a similar level of hypoxia, but for a much longer period, the entire 3-month life cycle, which may explain why they observed differential methylation in the sperm and we did not. If alterations to DNA methylation are not responsible for the changes in gene expression observed in our study, then perhaps histone modifications or miRNAs are at play. miRNA in medaka have also been identified as a potential underlying mechanism of transgenerational testis impairment induced by hypoxia by targeting genes associated with stress responses, cell cycle, epigenetic modifications, sugar metabolism, and cell motion [80].

Conclusion
Our study demonstrates that paternal exposure to hypoxia is able to endow offspring with a higher resistance to hypoxic conditions and that upregulation of genes at the hemoglobin gene enhancer locus on chromosome 3 may explain these phenotypic effects. This study is the first to demonstrate that paternal exposure alone can mediate these changes in phenotype and gene expression in response to hypoxia. These findings suggest that paternal inheritance may be providing offspring with valuable environmental specific information (an epigenetic memory [81]) about a hypoxic environment that could increase their survival, though it is unclear how this information is transmitted as we did not detect any differential methylation in the sperm of control and hypoxia treated males. Regardless, the variability in tolerance to acute hypoxia that we observed between hypoxia families suggests that progeny of different genotypes respond to hypoxia in different ways, phenotypically (i.e., an epigenotype x environment interaction). In other words, environmental specific information is transferred through some genotypes, but not others, suggesting that pre-acclimation to hypoxic conditions may be driven by an interaction between the genotype, the epigenotype and the environment.

Fish husbandry
Breeding and husbandry took place within the Otago Zebrafish Facility (OZF), a temperature-controlled facility maintained at 25-27 °C, pH 7-7.8, and conductivity 300-500 μS. Fish were maintained in a Tecniplast re-circulating system (Tecniplast, Varese, Italy) under a 14:10 light to dark photoperiodic cycle, with 30 min of simulated dawn and dusk at the start and end of each day. Acute hypoxia assays took place in the Zoology Department, where room temperature was controlled at 25 °C, with a 13.5-h (0700-2030 h) light cycle with 30 min of simulated dawn and dusk. Zebrafish were fed twice daily with dry food (ZM000-400, size-dependent) and once daily with live rotifer (Brachionus spp. days 5-10 post fertilization) or artemia (Artemia salina; days 10+ post fertilization). All animals were collected and maintained according to the standards of the Animal Ethics Committee for the University of Otago, New Zealand (protocol no. AEC 44/16).

Hypoxia exposure
In November 2016, 9-month-old male zebrafish (AB wild-type; n = 20/treatment) were exposed to hypoxic conditions (10.91-12.33 kPa pO 2 ) or control conditions (20-21 kPa), for 2 weeks, in line with Ho and Burrgren (2012). Two glass tanks (36Wx29Lx26.5H cm) were separated into three zones (12 × 29 cm), with 10 fish in each outer compartment (n = 20 per treatment) and two fine bubble diffusers and a filter positioned in the middle compartment. The oxygen concentration of the treatment tank was maintained by using an OxyGuard Mini probe (OxyGuard International, Denmark) and oxygen controller (Model PR5714, PR Electronics, Denmark) that were connected to nitrogen and air cylinders (BOC Gas Supplies, Food Fresh grade). An on/off relay output from the controller actuated a solenoid-controlled flow of compressed nitrogen or compressed air (BOC Gas Supplies, New Zealand), to maintain the system at 53.1-60% air saturation (= 10.91-12.33 kPa or 4.38-4.95 mg/L) level [59]. Normoxic conditions in the control tank (> 95% saturation) were maintained by continually passing air through the bubble diffusers, connected to an air pump. The oxygen concentration in the tanks was confirmed with a YSI 85 probe (YSI, Inc., Ohio, USA). Both tanks were siphoned for waste and received 10% water changes every 3 days. Calibration of the oxygen probe was checked daily and recalibrated as necessary. There was no mortality in either control or the hypoxia treatment.

F1 families, hypoxia tolerance assays, and analysis
Progeny from five families per treatment (n = 5 hypoxia, n = 5 control) were generated using a modified in vitro fertilization (IVF) protocol [65,82] 7 days after parental exposure finished (note that we generated IVF families on day 1 after exposure as well, but there was a technical issue with the sperm samples for downstream analyses). Briefly, eggs and sperm were collected using abdominal massage of females and males that were anesthetized by immersion in 0.01% Tricaine methanosulfonate-half the sperm (0.5-2 μL/male) was used for IVF, and the other half was stored for subsequent DNA extraction. One male and one female were used to create each F1 family.
Offspring were reared to 20-21 dpf when n = 8 offspring per family (except family C5, where n = 6; total n = 77) were challenged by acute hypoxia assays in small acrylic chambers (100 mm H × 50 mm × 50 mm). The 20-21 dpf was chosen as the larvae were then large enough to be reliably tracked with automated tracking software (see below), rather than timing manually by eye. Automated tracking increases repeatability and reduces observer bias. For the acute hypoxia assay, a 0-1 kPa oxygen level was achieved by continually passing compressed nitrogen through a bubble diffuser for at least 10 minutes before assays and continuously passing the nitrogen through a valve in the top of the chambers during the assay. The oxygen concentration fluctuated between 0 and 1%, monitored using a fiber-optic oxygen probe (Foxy OR-125) attached to an Ocean Optics ® USB 2000 spectrophotometer with USB-LS-450 light source and the manufacturer's software (OOI Sensors). Each fish was filmed for 240 s and recorded with a GoPro Hero 3+ camera. A total of 77 videos (n = 37 control, n = 40 hypoxia) were tracked using EthoVision XT behavioral tracking software, version 11.5 (Noldus Information Technology, Netherlands). Fish were then immediately euthanized by submersion in ice and stored in RNAlater (Invitrogen). One offspring from family C1 stayed at the bottom the entire assay, so this individual was removed from further analyses.
Resistance to acute hypoxia, defined as the first time that progeny lost equilibrium for 3 seconds or more during hypoxia assays [25], and the loss of equilibrium frequency were analyzed using R v. 3.5.1 [83]. Loss of equilibrium > 3 s was initially modeled using a Gaussian generalized linear model, and loss of equilibrium frequency as a Poisson generalized linear model, with treatment (control vs. hypoxia) as a fixed effect. To account for multiple fry per family being tested, we also ran the models incorporating treatment ID family as a random effect, using a Gaussian linear mixed effects model for time to loss of equilibrium > 3 s and a Poisson linear mixed effects model for the loss of equilibrium frequency. Significant heterogeneity was obvious in the time to loss of equilibrium, so a linear mixed model incorporating heteroscedasticity was compared to the model assuming homogeneity, by incorporating the variance in Treatment into the model [84]. Offspring from families H1-H3 and C1-C3 were chosen for further analyses-H1 and H2 were chosen as they showed the greatest phenotypic effect and all other families were chosen at random.

RNA-Seq library preparation and analysis
Total RNA from six whole 20-21 dpf offspring (3 control offspring, 3 treatment offspring) was extracted using a Zymo Duet extraction kit (Zymo, NZ). The integrity of RNA samples was determined using an Agilent RNA 6000 Nano chip on an Agilent 2100 Bioanalyzer to check that the samples had an RNA Integrity Number (RIN) value of 8-9. Total RNA concentration was measured by Qubit 2.0 Fluorometer (Qubit RNA HS Assay Kit, Life Technologies). Samples were sent to the Otago Genomics and Bioinformatics Facility at the University of Otago, under contract to New Zealand Genomics Limited, for library construction and RNA sequencing. Messenger RNA sequencing libraries were prepared using the Illumina TruSeq Stranded mRNA sample preparation kit (Illumina), as per the manufacturer's instructions. RNA sequencing was performed on the Illumina HiSeq2500 (Illumina, USA) machine with single-ended 100-bp reads generating 17.4-20.0 million reads per sample.
Low quality reads and remaining adapters were trimmed using Trim Galore v0.6.4 (https:// www. bioin forma tics. babra ham. ac. uk/ proje cts/ trim_ galore/) in a two-step process. First, sequencing adaptors were removed, 10 bp were hard-trimmed from the 5′ end to account for sequence bias produced by PBAT library preparation, and last, low-quality ends from reads (PHRED score < 20) were removed [85]. The reads were then aligned to the zebrafish genome (GCRz11) using HISAT2 v2.2.0 [86], informed with the GCRz11.99 annotation. Expression was summarized sample by sample at the gene level using featureCounts v2.0.0 [87]. The analysis of differential expression was conducted in R v.3.5.0 [83] using the DESeq2 package v1.24.0 [88] without the independent filtering option implemented in the "Results" function (see https:// doi. org/ 10. 5281/ zenodo. 68796 99). Differentially expressed genes between offspring of fathers exposed to hypoxia and controls were extracted after correcting for multiple testing using False Discovery Rate cut-off of q = 0.05 [89]. In order to test for over/under-representation of biological pathways which differentially expressed genes are involved, enrichment of Gene Ontology (GO) analysis terms was performed using Gorilla online accessed on August 4 2020 [90].

Whole genome bisulfite sequencing and analysis
DNA from six parental sperm samples (parental sperm of H1, H2, H4, C3, C4, C5 F1 families-note that only H1, H2, and C3 match the RNASeq samples, as there were some technical issues with library prep) was purified using a modified magnetic bead method [91]. WGBS was undertaken using an adapted modified post-bisulfite adaptor tagging (PBAT) method [67,92]. Briefly, bisulfite treatment was performed according to the EZ Methylation Direct Mag Prep kit (Zymo, D5044) instruction manual. Bisulfite treatment was performed before adaptor tagging, enabling simultaneous conversion of unmethylated cytosines, DNA fragmentation, and improving library preparation efficiency. Sequencing primers were added using random heptamer primers, and finally, sample-specific indexes and sequences required for Illumina flow cell binding were added by PCR. Library integrity was assessed by agarose gel electrophoresis and a fragment analyzer (Agilent) and sequenced on eight lanes (one flow cell) of an Illumina HiSeq using HiSeq2500 V4 sequencing of 100 bp single ended reads (in combination with 13 other samples, which were part of another study).
Raw reads were trimmed in Trim Galore v0.6.4 as previously mentioned. Read mapping was performed using Bismark v0.22.3 [93] with the option --pbat specified. Zebrafish genome version 11 (GRCz11) was used as reference. BAM files were deduplicated and reports containing methylation base calls were generated using the deduplicate_bismark and bismark_methylation_extractor scripts, respectively. The non-conversion rate during the bisulfite treatment was evaluated by calculating the proportion of non-CG methylation; by this measure, all libraries must have had a bisulfite conversion efficiency of at least 98.5% (Additional file 3: Table S3).
CG methylation calls were analyzed in SeqMonk v1.47.0 (www. bioin forma tics. babra ham. ac. uk/ proje cts/ seqmo nk). To analyze methylation at gene level, probes with a minimum of 5 methylation calls were generated and the percentage methylation measured as number of methylated calls/ total calls (mean CpG call depth ranged from 1.68-2.06 across the 6 samples). CpG islands (CGIs) were identified using Gardiner-Garden & Frommer's criteria and previously published datasets. For the former, 200-bp windows moving at 1-bp intervals were considered CGIs if the Obs/Exp value was greater than 0.6 and a GC content greater than 50%. For the latter, data obtained by biotinylated CxxC affinity purification (Bio-CAP) and massive parallel sequencing was used to identify non-methylated CpG islands [69]. Mean CpG depth at CpG islands ranged from 1.74-2.22 across the 6 samples.
Coupling between methylation and gene expression in differentially expressed genes (DEG) and at genome level was interrogated using methylation levels at transcription start sites (TSS). TSS were defined as 200 bp centered on the first nucleotide of an annotated mRNA, and a threshold of at least 20 methylation calls to be included in the analysis (mean TSS quantification ranged from 3.53 to 4.39 across the 6 samples). DEG were divided into underexpressed and overexpressed, whereas coupling at genome level was assessed dividing gene expression levels into quartiles. Custom annotation tracks were generated using Gviz v1.28.3 [94].
Finally, a conventional differential methylation analysis using the DMLtest function of the DSS package [95]. Briefly, DMR were identified by contrasting control and hypoxia samples using the DMLtest and callDMR functions in the DSS package, with the following parameters: smoothing = TRUE, smoothing.span = 500, minimum length DMR = 50 bps, minimum number of CpG sites = 5, delta DMR = 0 and p.threshold DMR = 1e−5. DMRs were annotated with the closest gene (threshold 2000 bp) using SeqMonk.